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SPATIAL MODEL OF CHEMOTHERAPY 



Abstract. In this paper we consider chemotherapy in a spatial model of tumor 
growth. The model, which is of reaction-diffusion type, takes into account 
the complex interactions between the tumor and surrounding stromal cells by 
including densities of endothelial cells and the extra-cellular matrix. When no 
treatment is applied the model reproduces the typical dynamics of early tumor 
growth. The initially avascular tumor reaches a diffusion limited size of the 
order of millimeters and initiates angiogenesis through the release of vascular 
endothelial growth factor (VEGF) secreted by hypoxic cells in the core of 
the tumor. This stimulates endothelial cells to migrate towards the tumor 
and establishes a nutrient supply sufficient for sustained invasion. To this 
model we apply cytostatic treatment in the form of a VEGF-inhibitor, which 
reduces the proliferation and chemotaxis of endothelial cells. This treatment 
has the capability to reduce tumor mass, but more importantly, we were able to 
determine that inhibition of endothelial cell proliferation is the more important 
of the two cellular functions targeted by the drug. Further, we considered the 
application of a cytotoxic drug that targets proliferating tumor cells. The 
drug was treated as a diffusible substance entering the tissue from the blood 
vessels. Our results show that depending on the characteristics of the drug it 
can either reduce the tumor mass significantly or in fact accelerate the growth 
rate of the tumor. This result seems to be due to complicated interplay between 
the stromal and tumor cell types and highlights the importance of considering 
chemotherapy in a spatial context. 

1. Introduction. Cancer is a complex multiscale disease and although advances 
have been made in several areas of research in recent years the dynamics of tumor 
growth and invasion are still poorly understood. It has been recognized that tackling 
this problem requires the collaborative effort of scientists from several disciplines 
such as genetics, cell biology, biological physics and mathematical biology to name 
but a few. This paper introduces a mathematical model of tumor invasion that was 
developed in precisely such a multi-disciplinary effort. Far from claiming that this 
type of collaboration is groundbreaking, we would still like to highlight the potential 
behind such a inter-disciplinary collaboration in cancer research. 

The second annual workshop of the Vanderbilt Integrative Cancer Biology Center 
(VICBC) took place from July 18-21, 2006 at Vanderbilt University in Nashville, 
Tennessee. During that workshop, our group of 16 biologists, mathematicians and 
biomedical engineers developed a mathematical model for tumor invasion. The 
model focuses on the complex interactions between different types of cells (nor- 
moxic, hypoxic, endothelial cells etc.) and their interaction with the surrounding 
extracellular matrix. The mathematical model is fully continuous with respect to 
time and space and consists of partial differential equations of reaction-diffusion 
type. The aim of this paper is to study cancer chemotherapy in a spatial context. 

The plan of the paper is as follows: First we discuss the biological background 
of the problem and discuss several aspects of tumor invasion and therapy includ- 
ing previous mathematical modelling. We then introduce our model and present 
numerical simulation results. We discuss our results in the context of other mathe- 
matical models and clinical applications and we outline possible directions for future 
research. An appendix describes the non-dimensionalization of the model. 

2. Biological background. 

2.1. Cell proliferation and survival. Cancer is regarded as a group of diseases 
characterized by abnormal cell growth and division, leading to tumors and eventu- 
ally metastasis. When the disease is not present, the healthy cells in the tissue have 
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the ability to proliferate and nutrients are kept at normal levels. They follow an 
ordinary cell cycle and upon growth and division create an increased cell density; 
when a certain threshold is met, the cells push on each other, creating a "crowd- 
ing effect." In this case, normal cells eventually become quiescent or may even go 
into apoptotis (programmed cell death). However, when normal cells become tu- 
morigenic, the cells will grow into three dimensional structures to overcome some 
of the crowding effect. The healthy balance between normal and apoptotic cells is 
disrupted and inside the tumor arises a new type of cell capable of sustaining itself 
with a lower nutrient level than the normal cells: the hypoxic cell. One of the most 
important roles played by hypoxic cells within a tumor is their ability to recruit 
new vasculature [16] to provide an additional supply of nutrients, such as oxygen, 
glucose, and growth factors. Oxygen is one of the most important nutrients in the 
study of cancer. It has been widely studied and documented with respect to its 
levels for various tumors due to the ease of measurement, and we have thus used it 
as our nutrient of choice for the work presented here. 

2.2. Unbiased and biased cell migration. One of the key aspects of cancer 
invasion is cell motility. Tumor cells can move within their environment by both 
undirected and directed mechanisms. Unbiased cell migration maintains no rec- 
ognizable pattern to the movement while cells undergoing biased migration follow 
a particular pattern, generally along a gradient, driven by physical or chemical 
means. This mathematical model will describe both types of movement using four 
mechanisms of motility. Diffusion, haptotaxis, crowding (as driven by tumor and 
non-tumor cells), and chemotaxis (as observed in angiogenesis). 

Diffusion is unbiased cell movement similar to the phenomenon of Brownian 
motion. This type of movement describes each cell's ability to randomly move 
within the surrounding environment. This random motility is applicable to all 
cell types included within our model. As the tumor becomes progressively more 
invasive, the tumor cells migrate toward dense areas of extracellular matrix (ECM) 
[9, 45]. This directed movement towards higher density of ECM proteins is known 
as haptotaxis. 

As the tumor cells proliferate, the tumor expands outward. As pressure increases 
from this proliferation, the cells are pushed into a new area [28, 35]. This pressure- 
driven movement is described as crowding. An example of this crowding effect 
involves colorectal cancers. Adenomatous tissue forms and invades within the villi 
of intestinal epithelium. As the space within the villi fill, the villus epithelium tears, 
and the adenomous tissue extends into the gut [13]. The tumor cells will continue 
to fill available space until crowding creates enough pressure to cause compression 
or damage to surrounding tissue; thus, crowding perpetuates the aggression of the 
tumor. 

2.3. Hypoxia-driven angiogenesis. As the cells proliferate and the tumor grows 
in size, nutrients and oxygen become the most critical factors for cell viability. 
When a small avascular tumor exceeds a certain size (with a critical diameter of 
approximately 2 mm), simple diffusion of oxygen to highly metabolizing tissues 
becomes inadequate, and hypoxia will be induced in the cell population lacking 
nutrients. This will, in turn, trigger the production and release of numerous tumor 
angiogenic factors (TAFs) to stimulate angiogenesis [30]. Pathological angiogenesis 
is hence a hallmark of cancer. Vascular endothelial growth factor (VEGF) is the 
most prominently induced molecule and plays a predominant role in angiogenesis 
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[29]. This process not only affects the blood supply to meet the metabolic demands 
of the growing mass of cells, but also represents the track to deliver anti-cancer 
drugs to all regions of a tumor in effective quantities [20] . 

2.4. Chemotherapy. Chemotherapy technically refers to the use of drugs that 
selectively target the cause of a disease, but is now almost synonymous with the use 
of drugs to treat cancer. The mechanisms of these drugs are very diverse, but they 
usually target cell proliferation by for example damaging the DNA or disrupting the 
mitotic spindles, which are essential for successful cell division. Chemotherapeutic 
drugs can be divided into three classes. Class I drugs are non-phase specific and 
affect both quiescent and proliferating cells; class II drugs are phase-specific and 
only target cells in a given phase; and finally class III drugs only kill cells that are 
proliferating. One example of a class III drug is Docetaxel [43], which induces cell 
death during mitosis by stabilizing micro-tubules and thus hindering the formation 
of functional mitotic spindles. In this paper we will focus on the modelling of class 
III drugs and investigate how the characteristics of the drug influences the efficacy 
of the treatment. 

The drugs mentioned so far are considered cytotoxic as they actively kill tumor 
cells. Another family of drugs we will consider are cytostatic drugs, which act by 
inhibiting cell division or some specific cell function in the cancer cells or host cells 
directly involved in tumor invasion. One example of this type of drug is Tamoxifen 
[41] used in the treatment of breast cancer, which binds to estrogen receptors on 
the cancer cells and therefore inhibits transcription of estrogen-responsive genes. 
Another example is Avastin [3], which does not target the cancer cells directly but 
instead affects the surrounding endothelial cells reducing the formation of new blood 
vessels supplying the tumor with nutrients. Finally, some drugs have been shown 
to have both cytostatic and cytotoxic effects, such as Lapatinib [37]. 

3. Mathematical models of tumor growth and treatment. Many researchers 
have proposed mathematical models for tumor growth and treatment. Let us men- 
tion here [8, 28, 7, 24, 58, 12, 34] among others. A recent book in the Springer 
Lecture Notes in Mathematics series [31] contains further interesting survey arti- 
cles. In these works the authors have investigated phenomena such as invasion, 
angiogenesis, and cytotoxic chemotherapy. 

Mathematical models have considered random motility as a major mechanism 
of cellular motion [48]. With continuous diffusion well-understood mathematically, 
this assumption forms the backbone of numerous models, including the majority of 
those that will be mentioned here. Random motility has been studied in the context 
of gliomas by Swanson [64, 36], who has done extensive comparisons with patient 
data. It has also been used to motivate the probabilistic movement of individual 
cells in models by Anderson and Chaplain [6] and Anderson [8]. Random motility 
is especially convenient because it is mathematically identical to the diffusion that 
governs chemical motion. This allows the inclusion of chemical species, such as nu- 
trients and tumor angiogenic factors, while remaining in the coherent mathematical 
framework of reaction-diffusion equations [38] . 

The model presented in this paper examines cell types in terms of spatially- 
dependent cell densities, and chemical species in terms of spatially-dependent con- 
centrations. This does not permit analysis of small-scale structures, particularly 
the vascular networks simulated in many angiogenesis models [6, 8, 22, 24, 40]. 
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However, it does allow for a large-scale perspective, and gives us a broader range of 
options in the interactions we consider. 

Tumor-induced angiogenesis has been the subject of considerable attention in 
mathematical models. Tumor angiogenic factor (TAF) plays a central role in driving 
such models, providing the means for a growing tumor to initiate the growth of new 
vasculature. The initial model of Chaplain and Stewart [25] uses a free-boundary 
problem to describe the interaction of TAF and endothelial cells implicitly. Since 
then, however, chemotaxis has become the normal way of modelling this interaction, 
with endothelial cells moving up the gradient of TAF concentration [26, 6, 8, 52]. 
A number of models also include haptotaxis of endothelial cells up a gradient of 
extracellular matrix [6, 22, 24, 52]. We have focused on chemotaxis as the dominant 
mechanisms of endothelial cell movement. 

Hypoxia is widely believed to be crucial to angiogenesis, with hypoxic cells se- 
creting TAF. A number of models have considered the transition of tumor cells 
between various states, such as Adam and Megalakis [2] and Michelson and Leith 
[46], and the diffusive model of Sherratt and Chaplain [60]. But in these models 
hypoxia is only included in the sense of quiescent, non-proliferating cells. On the 
other hand, vessel recruitment stimulated by TAF-secreting hypoxic cells has been 
modelled in a hybrid-discrete setting [5]. 

Chemotherapy has been subject to extensive mathematical modelling. These 
efforts have mainly focused on drug resistance [1, 27, 68], cell cycle specific drugs 
using age-structured populations [32, 37], and optimal treatment scheduling [10, 
49, 54]. Most previous models of chemotherapy have neglected the spatial effects 
and therefore disregarded an important component in tumor growth. Low levels of 
nutrient within the tumor can cause a large fraction of the tumor cell population 
to be quiescent and these cells are less responsive to chemotherapy. Because this 
is a spatial effect caused by the diffusion of nutrients from the blood vessels to the 
tumor, we believe that it is important to include space in a mathematical model 
of chemotherapy. This has been investigated in both a hybrid-discrete approach [5] 
and a continuous setting [50], but the model presented in this paper is different in 
that it takes into account the tumor- host interactions in a more detailed manner. 



4. The mathematical model. Our mathematical model is based on partial dif- 
ferential equations for the densities of all cells and chemical factors involved. Thus 
it is fully continuous with respect to time and space. Let t and x denote time and 
space, respectively. The complete set of dependent variables are listed in Table 1. 
The density of all types of cells and matrix combined is 

V = h J ra J rn J rf J rm. 

The equation for oxygen concentration is 

dw d^w 

— (x,t) = D w-q^2 + a wm(wm ax - w) - (3 w (n + h + m)w - >y w w, (1) 

where D w is the coefficient of oxygen diffusion, w max is the maximum oxygen den- 
sity, and (3 W is the uptake rate of oxygen by normoxic, hypoxic and endothelial cells. 
Oxygen is a diffusible substance that is provided at a rate a w by the vasculature, 
whose amount in the unit volume is assumed to be proportional to the endothelial 
cell density. This modelling approach is found for instance in [8]. The reason for 
the choice of the source term a w w(w rnax — w) is that, at high environmental levels 
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of oxygen, less oxygen is released through the vessel walls (i.e. by the vasculature). 
Finally oxygen decays at a rate j w . 
The normoxic cells are governed by 

dn , » d f /ry , ~ ^ .dn\ d ( df\ 

m {x > l) = & ( (D " max{n " v - 0} + Dm) d^)~&c \^ n i x ) (2) 

+ a n n{v max - v) - a h H(w h - w)n + —a h H(w - w h )h. 

Normoxic cells possess a "background" random motility D m . For simplicity, we 
assume that this random motility is constant, although recent work by Pennacchietti 
et al. [56] indicates that cells under hypoxic conditions can become more mobile. 
On top of this, a concentration of normoxic cells above a threshold v c adds to 
the dispersion of these cells through crowding- driven motion represented by the 
nonlinear diffusion term D n max{n — v c ,0}n x . Such a pressure driven motility 
was first proposed in a paper by Gurtin and MacCamy [35]. On the first line 
of equation (2) we find the term responsible for the haptotactic movement up a 
gradient of extracellular matrix. The haptotactic coefficient is denoted by \ n . The 
second line describes the gain and loss terms for normoxic cells. First, we assume 
a logistic growth with rate a n . The growth levels off in regions where the sum of 
all cells and matrix approaches the maximal density v max . In those regions where 
the concentration of oxygen drops below a certain critical value Wh, normoxic cells 
enter the hypoxic class at a rate ah- Here TL denotes the Heaviside function which 
is 1 for positive arguments and zero otherwise. The transition process is reversible, 
and hence there is an influx from the hypoxic class at a reduced rate j^^h (to our 
knowledge, this reduction factor is not available in the literature, we choose the 
value yq). We assume that it takes some time for a cell to resume the cell cycle 
after it has been recruited from a quiescent state. Note that the two processes are 
mutually exclusive. 

The equation for hypoxic cells is 

^(x, t) = a h TL{w h - w)n - ^a h H(w - w h )h - faH(w a - w)h. (3) 

The first and the second terms are dictated by conservation of mass and correspond 
to terms in equation (2). The third term describes the transition of hypoxic cells 
to apoptotic cells at rate fa as the level of oxygen falls below a second threshold, 
w a < ^h- Hypoxic cells are less active in general due to reduced availability of 
oxygen and other nutrients and we assume that lack of energy causes them to be 
immobile. 

The equation for apoptotic cells is given by 

0(1 

— (x,t)=p h H(w a -w)h. (4) 

The first term corresponds to the third term in equation (3). 
The equation for endothelial cells is 

Otti f dim dcj \ 

~m M = dx [ Dm ^ ~ mXm d~ x ) + «™ m ^— " ( 5 ) 

For simplicity and to reduce the number of free parameters we assume endothelial 
cells possess the same random motility D m as normoxic cells. Endothelial cells 
respond via chemotaxis to gradients of angiogenic factor g and require the presence 
of angiogenic factor for proliferation. Proliferation is capped by the total density of 
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cells. The proliferation constant for endothelial cells is a m . In a previous version 
of the model we also incorporated a haptotaxis term (in response to gradients of 
extracellular matrix), however, this made little difference to the simulation outcomes 
discussed here. We therefore decided, for the sake of simplicity, to omit haptotaxis 
from the endothelial cell equation. The proliferation of endothelial cells is governed 
by the same volume exclusion constraint as the proliferation of normoxic cells. 
The equation for the extracellular matrix is 

Of 

J-( x ,t) = -f3 f nf. (6) 

Tissue matrix is degraded by the tumor cells according to a mass-action law with 
rate constant (3f. 

The equation for angiogenic factor (VEGF) is 

^(x,t) = D g ^+a g h- f3 g mg. (7) 

Angiogenic factor moves by standard linear diffusion with coefficient D g . Angiogenic 
factor is produced by hypoxic cells alone at rate a g: taken up by endothelial cells 
with a mass-action coefficient of (3 g . Endothelial growth factor is a survival factor 
as well, so the endothelial cells consume VEGF at all times, not just when they 
proliferate. 

This system of equations will represent our baseline model, but in order to include 
therapy with a cytotoxic drug we will extend our model in section 6.3. In the next 
section we will non-dimensionalize the model and discuss values for the parameters. 



5. Parametrization. In order to simplify the analysis and simulations of the 
model we non-dimensionalize the model using the characteristic scales shown in 
Table 2. Although key parameters in our model, such as the oxygen consump- 
tion rates and diffusion constants, have been reported in the literature, the non- 
dimensionalization makes it easier to estimate parameters that have not been ex- 
perimentally determined by comparing the influence of different processes to each 
other. A helpful article in this respect is [59]. The detailed non-dimensionalization 
is carried out in the appendix. 

The oxygen diffusion constant was set to D w = 10 -5 cm 2 s _1 , as in [8]. The 
oxygen consumption of cancer cells has been experimentally determined to f3 w = 
6.25 x 10 -17 mol cell -1 s -1 [21], and for simplicity we assume that all cell types 
have this consumption rate. The oxygen production rate is difficult to estimate as 
it depends on vessel permeability and the oxygen concentration in the blood. We 
therefore make a non-dimensional estimate and set it to a w = 1, of the same order 
as the consumption rate. The oxygen decay rate is set to the non-dimensional value 
7 ™ = 0.025 [8]. 

The random cell motility parameter D m = 10 -9 cm 2 s -1 has been estimated [17] 
and for simplicity we assume that the pressure driven motility D n is of the same 
magnitude. The density at which this motility occurs is estimated to be 80% of 
the maximum density v max . Parametrizing the haptotactic coefficient is difficult 
as we consider a matrix consisting of a variety of cells and macro- molecules (MM), 
while the haptotactic movement of cancer cells is only sensitive to certain molecules 
bound within the matrix. From experimental data we know that the concentration 
of these molecules within the matrix is in the range 10 -8 — 10 -11 M [66]. The 
dimensional haptotactic parameter has been estimated to x ~ 2600 cm 2 s -1 M -1 
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[6], and using this value we can at least give a non-dimensional estimate to the 
haptotaxis coefficient % n = 1.4 x 10 -4 [8]. 

By rescaling time by the average doubling-time of the cancer cells the growth 
rate is by definition a n = log 2. The conversion rate from normoxic to hypoxic cells 
at low oxygen pressure has to our knowledge not been experimentally determined, 
but a reasonable guess is that this occurs on the time-scale of hours. The rate is 
therefore estimated at a h = 2.8 x 1(T 5 s~\ which means that in hypoxic conditions 
half of the cells will turn hypoxic within log 2/2.8 x 10 -5 « 7 hours. The process 
of apoptosis is on the time-scale of minutes [67], but the time it takes before the 
cells commits to this decision is much longer and on the order of days [15]. The 
conversion rate from hypoxic to apoptotic cells is therefore estimated to be consid- 
erably smaller and set to flh = 5.6 x 10 -6 s _1 . The exact oxygen concentrations 
at which the above conversions occur is difficult to measure as it depends not only 
on the cell type, but also on other environmental factors such as acidity and the 
concentration of various growth factors. Instead these parameters can be estimated 
from oxygen concentrations measured in the necrotic centers of real tumors. Such 
measurements show that the oxygen concentration is approximately 0.5-30 % of 
that in the surrounding tissue [18]. From this we estimate the hypoxic threshold 
Wh to be 5 % of the background oxygen concentration and the apoptotic threshold 
w a to be slightly lower at 3 %. 

The chemotactic coefficient of the endothelial cells was set to \m — 2600 cm 2 s _1 
M _1 [6]. The growth rate of the endothelial cells was estimated non-dimensionally 
by comparing it to the growth rate of the cancer cells. Endothelial cells are known 
to have much longer doubling times than cancer cells [55], and we therefore set 
otm = & n /10 ~ 0.07. The matrix degradation rate (3f has not yet been measured 
experimentally and we therefore resort to an estimate of this parameter. If we 
assume that a cancer cell degrades a volume of matrix which is of the same order 
of magnitude as the size of the cell itself during one cell cycle, then the degradation 
rate can be estimated to /3f = 10 -13 cm 3 cell -1 s -1 . 

The diffusion coefficient of VEGF was set to D g = 2.9 x 10" 7 cm 2 s" 1 [6]. The 
production rate of VEGF by hypoxic cells was set to a g = 1.7 x 10 -22 mol cell -1 
s _1 , similar to the value used in [39] (0.08 pg cell -1 day -1 ). The uptake of VEGF by 
endothelial cells has not been experimentally determined and we therefore estimate 
it to be of the same order as the production rate and set /3 g = a g = 10. All 
parameters in the model are summarized in Table 3 in their non-dimensional and 
dimensional form, where appropriate. 

6. Results. We have implemented our model (13)- (19) using matlab (version 7.1, 
The MathWorks, Inc., Natick, MA). We used the function pdepe, which discretizes 
the equations in space to obtain a system of ordinary differential equations 1 . The 
system of equations (13)-(19) is solved in the domain Q = [0, 1] x (0,r ma J, where 
T~max will vary depending on the scenario we want to model. 
The model is supplied with the following initial conditions: 



n(&0) 
m(£,0) 



1.0, £e[0,l], 
0.93exp(-200£ 2 ), 
0.01, 

1 - ra(£, 0) - m(£, 0) - 0.05, 



(8) 



The MATLAB codes will be available from the corresponding author upon request. 
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all other initial data being zero. We consider an initial avascular tumor located at 
£ = in a tissue with a uniform background distribution of endothelial cells (blood 
vessels) and extracellular matrix. The oxygen concentration is initialized at its 
saturation value in the entire domain. We assume that all cells and chemical species 
remain within the domain, and therefore no-flux boundary conditions are imposed 
on dfl for all equations. This is clearly a simplification of the actual situation, but 
it will already give instructive results (see [23] for a similar situation) . First we will 
present simulation results from the baseline version of the model and in the following 
sections we will modify the model to incorporate several modes of treatment. 

6.1. Baseline. This case represents the baseline scenario when tumor growth only 
is limited by intrinsic constraints such as oxygen supply and the surrounding stromal 
tissue. The simulation lasted T max = 300 time steps (200 days) and the results of 
the simulation can be seen in Figure 1, which shows the spatial distribution of the 
normoxic cells (n), hypoxic cells (ft), endothelial cells (m) and the extracellular 
matrix (/) for four different time points of the simulation. Note that the density 
of the endothelial cells has been scaled up by a factor of 10 in all plots for better 
visualization. 

At the beginning of the simulation we see the initial tumor growing and invading 
the surrounding tissue by degrading the extracellular matrix, but when the tumor 
reaches a critical size (approx. 3 mm or after 250 division cycles) the oxygen supplied 
by the endothelial cells becomes insufficient and we observe the emergence of hypoxic 
cells in the core of the tumor (Figure 1). The hypoxic cells secrete VEGF, which 
triggers endothelial proliferation and chemotaxis towards the hypoxic region. This 
consequently increases the oxygen production and leads to sustained tumor invasion. 
The end result is a traveling wave like scenario with a leading edge of normoxic cells 
followed by hypoxic cells and finally by a bulk of apoptotic cells (not shown). The 
density of endothelial cells has a peak at the edge of the tumor, is zero inside the 
tumor and remains approximately at the initial value in the surrounding tissue. 
The density of endothelial cells at the tumor-host interface is five-fold higher than 
the initial value, which implies that an angiogenic response has occurred due to the 
hypoxia experienced by the cancer cells. 

The above scenario is firmly established as the typical scenario for early tumor 
growth [63, 57] and has also been reproduced in previous mathematical models 
[51, 5], but as many of the parameters by necessity were estimated rather than ex- 
perimentally determined we also investigated variations of these parameters within 
a reasonable range (data now shown). For example, increasing the oxygen con- 
sumption rate of the cells /3 W or the decay rate of oxygen j w leads to an earlier 
appearance of hypoxia. Increasing the reproduction rate of endothelial cells gives 
rise to a stronger angiogenic response. A decrease of the diffusion constant of the 
angiogenic factor D g results in a more persistent gradient of angiogenic factor and 
a wider spread of endothelial cells, without, however, a greater number of normoxic 
tumor cells. This outcome is also seen if the chemotactic sensitivity Xm of endothe- 
lial cells is increased. Future work will include a more systematic exploration of the 
parameter space and will indicate to experimentalists which parameters influence 
the tumor behavior most. 

6.2. Anti-angiogenic therapy. Ant i- angiogenic therapies are considered to be 
cytostatic in the sense that the drugs used are not toxic to the cells, but instead 
inhibit some mechanism essential for cell division or a specific function. Probably 
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the most well-known anti-angiogenic drug is Avastin (Bevacizumab) [3], which acts 
on endothelial cells by inhibiting the function of VEGF. The drug thus influences 
the endothelial cells' ability to react to the VEGF produced by hypoxic cells in the 
tumor and consequently is expected to reduce the angiogenic response and therefore 
inhibit further tumor invasion. 

We will model the effect of this type of drug by altering the behavior of the en- 
dothelial cells with respect to VEGF. The growth factor increases the proliferation 
and induces chemotactic movement of the endothelial cells, and we will model the 
effect of the drug by altering the parameters that correspond to these two processes. 
Avastin is usually administered intra-venously and therefore arrives at the tumor 
site through the blood vessels. As the endothelial cells in our model represent blood 
vessels these cells will be instantly affected by the drug, and there is therefore no 
need to introduce a new chemical species for the drug. Instead we can alter the val- 
ues of the parameters under consideration directly depending on the concentration 
of the drug in the blood. Avastin is usually administered in 1-2 week intervals, but 
because the half-life of the drug is very long (1-2 weeks) [11] the concentration of the 
drug in the blood can, as a first approximation, be considered constant. This im- 
plies that we can alter the parameter values at the start of the treatment and keep 
them constant during the entire treatment period. As mentioned before Avastin 
decreases the endothelial cell proliferation a m and the chemotaxis coefficient Xm> 
we therefore model the effect of the drug by reducing these two parameter by a 
factor 10 during the treatment, i.e. when r s < r < r e we let 

otm «m/10, 

Xm ~> Xm/10, 

where r s and r e are the start and end times of the treatment. In order to further 
investigate the effect of the drug we will decouple the inhibitory effect on VEGF and 
separately consider the impact of inhibiting the proliferation and the chemotaxis. 

The approach we have used is of course a highly simplified way to model the 
treatment. A more complete model would include the pharmacokinetics and phar- 
macodynamics of the drug and also treat the effect of the drug on endothelial cells in 
more detail, but our approach serves as a first attempt at incorporating the effects 
of an anti-angiogenic drug into a detailed spatial model of tumor-host interaction. 

The results of the anti-angiogenic drug on the system can been seen in Figure 
2, which shows the spatial distribution of normoxic, hypoxic and endothelial cells 
at (a) r = 250 and (b) r = 300. The drug was administered in the time interval 
200 < r < 300 (approx. 10 weeks) and in this simulation both effects of the drug 
were considered (i.e. both proliferation and chemotaxis). Up to the point where 
the treatment begins the dynamics are identical to that of the baseline scenario and 
are therefore not shown. When the treatment begins the endothelial cells become 
unresponsive to the VEGF and this has an obvious effect on the growth dynamics. 
The endothelial cells are now much more dispersed and not centered anymore at 
the leading edge of the hypoxic cell population (Figure 2a). This effect is even more 
pronounced at the end of the treatment (Figure 2b), when the endothelial cells are 
almost evenly distributed in the whole domain. This consequently decreases the 
oxygen supply to the normoxic tumor cells and leads to a lower cell density in the 
proliferating rim. The treatment could therefore be considered successful, but it 
should be noted that when the treatment ends the endothelial cells will regain their 
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proliferative and chemotactic abilities and this could potentially lead to re-growth 
of the tumor. 

Figure 4 shows the result of the treatment if the anti-angiogenic drug only affects 
the endothelial cell chemotactic sensitivity Xm (the duration of the treatment is as 
above). In this case the growth dynamics are quite similar to the baseline scenario, 
but on the other hand if only proliferation is affected (Figure 3) the dynamics 
are similar to the full treatment. This suggests that inhibiting the proliferation of 
endothelial cells is the more important of the two effects the drug has. 



6.3. Modelling cytotoxic therapy. We now expand our model to include therapy 
with a cytotoxic drug. The effects of this type of treatment has been investigated us- 
ing several mathematical models [53, 4, 1, 44, 68, 42, 65]. These studies have mostly 
focused on the efficacy of specific drugs, but we will instead model the effects of 
a general cytotoxic drug that affects cells in the proliferative state (i.e. normoxic 
cells), and investigate how the drug characteristics influence the outcome of the 
drug therapy. Many drugs used in cytotoxic chemotherapy, for example spindle 
poisons (paclitaxel) and drugs that interfere with the DNA synthesis (doxorubicin, 
fluorouracil) affect only proliferating cells. The concentration of the drug c is intro- 
duced as a new variable and it is described by a reaction-diffusion equation similar 
to that of oxygen, 

Oc O^c 

— (£,r) = D c~g^ + a c (r)ra(l - c) - 7 c c - k^ n nc. (9) 

Notice that this equation is already in non-dimensional form and that < c < 1. 
The drug is assumed to be delivered through blood infusion and therefore enters 
the tissue from the blood stream. The production rate of the drug is therefore 
proportional to the density of endothelial cells, which implies that regions with 
higher vascular density will experience a higher concentration of the drug. Instead of 
looking at a constant drug supply we include a scheduling of the drug by letting the 
production rate a c (r) be a time-dependent function, where a c (r) > corresponds 
to the drug being delivered to the tissue. This is of course a highly simplified 
picture of the true dynamics of drug delivery. A more realistic approach would be 
to also include the pharmacokinetics of the drug into the model [14], but as we will 
model a general cytotoxic drug we will use this simplified form and only consider 
the dynamics of the drug at the tumor site. The drug diffuses (for simplicity we 
assume the same diffusion constant D c as for oxygen) in the tissue and decays at a 
constant rate j c . It affects the normoxic cells alone and is consumed at a rate kj n 
when it kills normoxic cells. The normoxic cells are driven into apoptosis at a rate 
proportional to the drug concentration 

On 

g^(^ T ) = 7n^C, (10) 

with rate constant j n > 0. The dimensionless parameter k corresponds to the 
amount of drug needed to kill one unit volume of proliferating cells. For sake of 
simplicity, we have worked with k = 1. Finally, to balance the loss of normoxic 
cells, the equation for apoptotic cells is amended according to 

Ocl 

<^(£> r ) = ■■■+lnnc. (11) 
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The treatment schedule was fixed at 



5 



a, 



: c (r) = 100 exp(-4(r - (200 + 2k)) 2 ), 



which corresponds to a treatment that starts at r = 200 and consists of 6 infusions of 
the drug. Each infusion lasts approximately 1 cell cycle, and the interval between 
two infusion is 2 cell cycles. This treatment schedule might not be realistic for 
some drugs (e.g. due to toxicity), but in order to separate the effects of scheduling 
and drug characteristics we will use this basic schedule in our investigation. An 
application of this treatment schedule with parameters j n = 10 and j c = 0.1 can 
be seen in Figure 5, which shows the post-treatment density of normoxic, hypoxic 
and endothelial cells at r = 250. In this case the treatment is successful and the 
tumor mass is significantly reduced compared to the untreated baseline scenario 
(Figure 1). From this plot we can also observe that the reduction in normoxic cells 
is largest where the endothelial cells are located and that the surviving fraction of 
cancer cells now reside in a region of low vascular density. For other values of the 
drug parameters, in particular a reduced kill rate 7 n , the treatment is less successful, 
and it can even lead to enhanced tumor growth, as can be seen in Figure 6. This 
plot shows the post-treatment density of cells (at r = 250), and comparing it to 
the baseline simulation we can observe that the total number of normoxic cells is 
actually higher although chemotherapy has been applied. This is further illustrated 
in Figure 7, which shows the time evolution of the total tumor mass (normoxic 
plus hypoxic) for the baseline scenario and three different drug characteristics. For 
r < 200 the dynamics are identical as the initial conditions and governing equations 
are identical, but when the therapy is applied at r = 200 they diverge. In the initial 
growth stage (r < 150) we observe an approximately linear increase in tumor mass. 
This is followed by a sharp transition to a slowly decreasing tumor mass which is 
due to the build up of apoptotic cells in the center of the tumor (if the apoptotic 
cells are included in the tumor mass the linear growth rate is preserved). When the 
treatment is applied the dynamics clearly depend on the drug parameters and we 
observe that the treatment can decrease, but also increase the total tumor mass. 

In order to investigate this further and to fully characterize the impact of the 
kill and decay rate on the efficacy of the treatment we performed systematic mea- 
surements of the post-treatment tumor mass in the range [1,100] x [1,100] of j c 
and 7 n . This was done by running the simulation for 400 evenly distributed points 
in the parameter range and at r s = 300 measure the total mass of the tumor (nor- 
moxic plus hypoxic cells). The tumor mass was then normalized with respect to 
the untreated case. The normalized post-treatment mass is thus given by, 



where Mo is the untreated tumor mass. 

The results can be seen in Figure 8 and show that in a significant part of the 
parameter space the treatment actually increases the tumor mass. We can observe 
that the treatment is most successful for high kill rates and low decay rates, where 
the post-treatment tumor mass M « 0.2. On the other hand for intermediate kill 
and decay rates the tumor mass is higher than in the untreated case reaching a post- 
treatment mass of M w 1.5. For low kill rate and high decay rates we have M ~ 
1, which is expected as the untreated case corresponds to 7 n =0. The parameter 
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space can therefore be divided into three distinct regions: (i) Successful treatment 
(M < 0.75), (ii) Unaffected growth, and (iii) Accelerated growth (M > 1.2). This 
is illustrated in Figure 9, which shows that a large region in the parameter space in 
fact gives rise to accelerated growth. 



7. Discussion. The model presented here is similar to previous reaction-diffusion 
models of tumor invasion, but also contains a number of novel features that are 
worth mentioning. Firstly we have included density driven motion of the cancer 
cells, which acts by pushing cells out from regions of high density. This mechanism 
is potentially important as tumor tissue is known to be highly irregular and to have 
higher cell density than normal tissue. Although we have not seen any obvious 
effect of this mechanism, this might be different when considering the system in 
two dimensions, where this type of density driven motion is known to create front 
instabilities [47], which can cause fingering protrusions. Our approach of using a 
degenerate diffusion term to account for cell crowding is simpler and more transpar- 
ent than a volume fraction formulation with a force balancing term. This approach 
also meshes more easily with standard representations of haptotaxis and random 
motility. 

Secondly we have taken the extracellular matrix into account when considering 
the maximum density of cells. This implies that the cancer cells cannot move into 
regions where the matrix density is too high and thus imposes a growth constraint 
which clearly is present in real tissue. Further it implies that degrading the matrix 
does not only induce haptotaxis, through the establishment of matrix gradients, 
but also pure random movement, as the degradation frees up space that the cancer 
cells can move into. 

In the basic setting of our model, which represents unconstrained tumor growth 
and invasion, the simulations produce results that are in qualitative agreement with 
the typical scenario of early tumor growth: the tumor starts growing as an avascular 
tumor that reaches a critical diffusion limited size of a couple of millimeters. At 
this stage the tumor has outgrown its nutrient supply and cells in the center of the 
tumor become hypoxic due to the low oxygen concentrations. This triggers them to 
secrete VEGF which diffuses in the tissue and stimulates endothelial cells to form 
blood vessels towards the hypoxic site through proliferation and chemotaxis. As a 
result of this angiogenic response the tumor can continue to grow and invade the 
surrounding tissue. 

Although this does not give any new insight into the dynamics of tumor invasion it 
nevertheless gives us a reference with which we can compare the different treatment 
strategies we apply. This is important as it gives us the possibility to examine 
different types of treatment within a fully developed framework of tumor growth. 

The results from the simulations of ant i- angiogenic treatment show that it has 
the capability to reduce tumor mass and slow down tumor invasion. This is in 
qualitative agreement with experimental results [33], but the strength of modelling 
this treatment is that we have the capability to decouple the effects of the drug on 
endothelial cells. In doing so we observe that proliferation of endothelial cells seems 
to be more important in establishing sufficient angiogenesis than chemotaxis of en- 
dothelial cells. Reducing the chemotaxis alone seems to have little or no therapeutic 
effect, while reducing the proliferation of the endothelial cells has a similar effect 
as the full treatment. This is an interesting observation and suggests that drugs 
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that specifically target proliferation of the endothelial cells could be more efficient 
as ant i- angiogenic agents. 

In the case of the cytotoxic treatment we also made observations that are sur- 
prising and even counterintuitive. The fact that the cytotoxic treatment leads to 
accelerated tumor growth is quite unexpected, as even a drug with a low kill rate 
will reduce the number of normoxic cells in the tumor, and from a naive view-point 
this would reduce the tumor mass. The fact that we observe the opposite, an in- 
crease in the tumor mass, reveals that the dynamics of the growth are more complex 
than one would expect. The results suggests that there exists a complicated inter- 
play between the different cell types occupying the tumor and that disrupting the 
balance between these can have unexpected consequences. 

When the kill rate of the drug is low (or the decay rate is high) only a small 
fraction of the normoxic cells at the boundary of the tumor are driven into apoptosis. 
As the reduction in number of normoxic cells, at the tumor boundary, lowers the 
total oxygen consumption this leads to a conversion of previously hypoxic into 
normoxic cells. The interesting fact is that the number of normoxic cells after 
treatment is actually higher than before, which suggests that the local carrying 
capacity is increased in the accelerated growth regime. This can only occur if 
the local density of endothelial cells is increased, which implies that the treatment 
indirectly influences the vascular density at the tumor- host interface. This also 
highlights that cytotoxic therapy can have indirect effects on the micro-environment 
of the tumor and can even make the growth conditions more beneficial for the tumor. 
This observation might also shed new light on the efficacy of existing cytotoxic drugs 
and possibly guide the development of new drugs. It should be noted that a model 
used for this specific purpose would need to be properly experimentally parametrized 
and validated, however, at least the model presented in this paper is a first attempt 
at modelling the influence of tumor-host interactions on chemotherapy. The fact 
that the models takes into account the interactions between the tumour and the 
stromal cells means that the results need to be compared with in vivo rather than 
in vitro studies. One possibility is to compare tumor growth rates during cancer 
therapy with different drugs under the same delivery schedules. The decay/kill rates 
of these drugs could also be measured, which would make a comparison with the 
model results possible. 

The model originally contained further variables which we later disregarded for 
the sake of simplicity. In future extensions of this work we will consider different 
types of extracellular matrix (namely tissue matrix and clotting matrix) as well 
as reactive cells. Reactive cells are attracted by gradients of clotting matrix and 
can turn clotting matrix into tissue matrix. Eventually one may observe an encap- 
sulation of the tumor by fibrous tissue. Another simplification that was made in 
the present model is that matrix degrading enzymes are located in the cell mem- 
branes of normoxic cells alone. In future work we plan to include diffusible matrix 
degrading enzymes as a new dependent variable. Finally, we plan to take the nu- 
merical simulations to two-dimensional domains to allow for greater morphological 
complexity. 
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8. Conclusion. We have seen in the present paper that mathematical models with 
a spatial component allow for complex and sometimes counterintuitive growth dy- 
namics of tumors. Therefore, mathematical models serve as an ideal tool to evalu- 
ate the efficacy of cancer drugs. An interface between experimental biologists and 
mathematicians can streamline this process considerably. 

We built our model based on few basic assumptions summarized as follows: (i) 
initial tumor growth is only constrained by the surrounding tissue, (ii) the tumor 
mass rapidly outgrows nutrient supply, (iii) subsequent tumor cell death stimu- 
lates stromal response, (iv) new stroma develops vasculature and (v) tumor growth 
resumes. These assumptions are well supported by available literature [63, 57]. 
However, they are only a partial set of variables that one may consider in building 
a mathematical model of tumor growth. They were chosen because in our view 
they are a minimal set that can realistically capture the effects of cytotoxic and 
cytostatic drugs. Naturally, empirically determined parameters would have to be 
introduced in order to make the model entirely realistic. Nonetheless, there are 
novel features as well as conclusions of the model that are worth reporting at this 
stage. 

One novel feature of this model is the way tumor cell migration is modeled, i.e., 
tumor cells are driven by crowding, in addition to the more conventional unbiased 
diffusion and haptotaxis. This introduces a measure of realism because it is generally 
recognized that tumor invasion is partly due to a pressure build-up that pushes 
cells into the surrounding tissue. Further we have also taken into account the 
extracellular matrix when considering the maximum cell density in the tissue, which 
implies that the stroma restrains the growth of the tumor. 

Perhaps the most compelling novel feature of our model is chemotherapeutic 
drug delivery by a spatial source, i.e., endothelial cells. Cancer chemotherapy has 
been the subject of countless modelling efforts [1, 27, 68, 32, 37, 10, 49, 58, 54]. 
However, spatial aspects of chemotherapeutic drug delivery seem to have received 
less attention so far. This is probably because the main objective is to evaluate 
the ability of a drug to debulk a tumor by inhibiting growth kinetics or killing 
cells. In contrast, in our model we consider the effects of removing cells (by drug 
killing or growth inhibition) in specific locations of a growing tumor. That is, we 
view chemotherapeutic drugs as agents that have a localized spatial effect on the 
growth of the tumor. This means that the endothelial cells have a dual role in our 
model. They supply the tumor with essential nutrients, but also when therapy is 
applied, with detrimental drugs. The angiogenic response therefore becomes crucial 
as it is necessary for continued growth, but also is utilized for drug delivery, and 
again this highlights the importance of taking into account the complex dynamics 
of tumor-host interactions when modelling drug therapy. 
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Appendix. The length is re-scaled with the typical size of an early stage tumor 
L = 1 cm and time with T = 16 h the typical time of a cell cycle for cancer 
cells [19]. The cell densities are re-scaled with a maximum cell density of v max = 
10 -8 cell/cm 3 , which corresponds to each cell having a diameter of approximately 
20 /im [21]. The oxygen concentration is re-scaled by the background concentration 
Wmax — 6.7 x 10 -6 mol/cm 3 [8] and the VEGF concentration by the maximal 
concentration go = 10 -13 mol/cm 3 [6]. 
We scale the independent variables 

. x t 

The dependent variables are scaled as follows 

w n ~ h „ a 

-, h = , a — , 



9 tot 



Notice that the concentration of extracellular matrix / has the unit cells/cm 3 . To 
convert this into a concentration of the unit moles /cm 3 one has to multiply with 
a conversion factor of the unit moles/cell. This conversion factor can depend on 
the type of matrix molecules and the type of cells considered. Equation (1) for the 
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oxygen becomes 

dw . , TD W d 2 w _ „ . . _ P w Tv max ~ . „ „ 

^-(^ r ) ^2~^72" +^ Tv max^(l ~ w) + ft + m)W ~ J w Tw. 

OT L Ot, Wmax 

We define new dimensionless parameters 

j\ TD W ~ rp a /3wTv max „ 

Upon dropping the tildes everywhere for notational convenience, our new dimen- 
sionless equation becomes 

dw d 2 w 

— (f,r) D w~Q£2 + a wm(l - w) - f3 w (n + h + m)w - >y w w. 

A similar procedure was applied to the remaining equations (2)- (7) and by drop- 
ping the tildes they return to their original form with the only difference that the 
parameters now are in their non-dimensional form. The complete simplified non- 
dimensional system that we will solve numerically in the following sections is as 
follows, 

dw d 2 w 

— D ™~q^2 + ot w m{l - w) - p w (n + ft + m)w - >y w w, (13) 

m {x > l) = T X \} Dn max{n " v - 0} + D ^)-d- x [ Xnn ir x ) (14) 

+ a n n(v max - v) - a h H(w h - w)n + ^a h H(w - w h )h, 
dh 1 

— Or, t) a h H(w h - w)n - —a h H{w - w h )h - (3 h H(w a - w)h, (15) 

^(x, t) PhH(w a - w)h, (16) 
dm , d / dm dg\ , x 

df 

-±(x,t) = -p f nf, (18) 



\ ^ d 2 g 



(x,t) = D g ^-^ + o^ft - (3 g mg. (19) 
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t = 50 





T = 100 
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t = 300 
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Figure 1. The spatial distribution of normoxic, hypoxic, endothe- 
lial cells and extracellular matrix at r = 50, 100, 200, 300 for the 
baseline scenario. Notice that the density of endothelial cells is 
multiplied by 10 for better visibility. 
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Figure 2. The spatial distribution of normoxic, hypoxic and en- 
dothelial cells during and after anti-angiogenic therapy. The pro- 
liferation a m and the chemotaxis coefficient Xm of the endothelial 
cells are reduced by a factor of 10 while 200 < r < 300 (all other 
parameters are as in the baseline scenario). 
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Figure 3. The spatial distribution of normoxic, hypoxic and en- 
dothelial cells during and after ant i- angiogenic therapy when the 
proliferation a m of the endothelial cells is reduced by a factor of 10 
while 200 < r < 300 (all other parameters are as in the baseline 
scenario) . 
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Figure 4. The spatial distribution of normoxic, hypoxic and en- 
dothelial cells during and after ant i- angiogenic therapy (at times 
r = 250 and r = 300) when the chemotaxis coefficient Xm of the 
endothelial cells is reduced by a factor of 10 while 200 < r < 300 
(all other parameters are as in the baseline scenario). 



t = 250 t = 300 




Figure 5. The spatial distribution of normoxic, hypoxic and en- 
dothelial cells after chemotherapy has been applied with the drug 
parameters 7 n = 10 and j c = 0.1 at times r = 250 and r = 300. 
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Figure 6. The spatial distribution of normoxic, hypoxic and en- 
dothelial cells after chemotherapy has been applied with the drug 
parameters j n = 1 and j c = 1 at times r = 250 and r = 300. 
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Figure 7. The time evolution of the total tumor mass (normoxic 
plus hypoxic cells) for the baseline scenario and three different 
chemotherapy treatments. It can be observed that chemotherapy 
not only has the capability to decrease, but also to increase the 
growth rate of the tumor. The decrease in tumor mass prior to 
treatment is due to the emergence of apoptotic cells in the center 
of the tumor and the linear growth rate would be preserved if the 
apoptotic cells were included in the total tumor mass. 



SPATIAL MODEL OF CHEMOTHERAPY 



23 




Figure 8. Post-treatment normalized tumor mass as a function 
of the kill rate (jy n ) and decay rate (7 C ). The untreated case corre- 
sponds to M = 1. 
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Figure 9. The division of drug parameter space into three dis- 
tinct regions: (i) successful treatment (ii) normal growth and (hi) 
accelerated growth. Here successful treatment is defined by a nor- 
malized tumor mass M < 0.75 and accelerated growth by a nor- 
malized tumor mass M > 1.2 (see equation (12) for the definition 
of M). 
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Notation 


Concentration of . . . 


w 


nutrient 


n 


normoxic cells 


h 


hypoxic cells 


a 


apoptotic cells 


m 


endothelial cells 


f 


extracellular matrix 


9 


VEGF 


c 


cytotoxic drug 



Table 2. The scaling factors for the nondimensionalization procedure. 



Parameter 


Value 


Meaning 


T 


16 h 


cell cycle time 


L 


1 cm 


typical length scale 


Vmax 


10 8 cells/cm 3 


maximum density of cells 


Wmax 


6.7 x 10" 6 moles /cm 3 


saturation level of oxygen 


9tot 


10 — 13 moles /cm 3 


maximal concentration of VEGF 



Table 3. Model parameter values in both their nondimensional 
and dimensional forms. 



Parameter 


ND-value 


D-value 


Reference 


D w 


0.58 


10" 5 cm 2 s" 1 


[8] 




0.57 


6.25 x 10- 17 mol cell" 1 s' 1 


[21] 




1 








0.025 




[8] 




5.8 x 10- 5 


10~ 9 cm 2 s- 1 


[17] 


D n 


5.8 x 10- 5 




[17] 


Vc 


0.8 






Xn 


1.4 x 10~ 4 




[8] 


U n 


log 2 






Cth 


1.6 


2.8 x 10~ 5 s- 1 




Ph 


0.32 


5.6 x 10~ 6 s- 1 


[15] 


Wh 


0.05 


3.4 x 10 _T mol cm -3 


[18] 


Wa 


0.03 


2.0 x 10 _T mol cm -3 


[18] 


Xm 


2.1 x 10~ 6 


2.6 x 10 3 cm 2 s -1 M _1 


[62] 




0.07 




[55] 


Pf 


0.5 


lO -13 cm 3 ceir 1 s- 1 






0.02 


2.9 x 10~ 7 cm 2 s -1 


[61] 


OLg 


10 


1.7 x lO" 22 mol cell" 1 s" 1 


[39] 


P 9 


10 







